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ABSTRACT 


This  report  is  concerned  with  the  development  of  a  fast  computer 
method  for  evaluating  the  albedo  integral.  This  integral  defines  the 
illumination  on  an  arbitrarily  oriented  surface  element  at  any  point 
in  space  about  a  diffusely  reflecting  sphere.  It  enters  the  calculation 
of  simulation  control  parameters  in  the  Arnold  Engineering  Develop¬ 
ment  Center  Aerospace  Environmental  Chamber  (Mark  I).  The  semi- 
numerical  method  developed  here  is  faster  than  ordinary  numerical 
integration  by  a  factor  of  about  ten.  A  typical  computer  program, 
which  formerly  required  about  thirty  minutes,  now  produces  the  same 
results  in  under  four  minutes. 
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NOMENCLATURE 

Albedo,  fraction  of  solar  radiation  reflected  by  earth 

Horizon  circle 

Terminator  curve 

"Target  plane  cut"  curve 

Surface  element  on  albedo  source 

Target  surface  element 

Altitude 

Solar  constant,  intensity 

Intensity  of  target  illumination 

Orbit  angular  momentum  vector 

Vector  normal  to  target  surface 

Solar  node  vector 

Target  position  vector 

Radius  of  the  albedo  source  sphere 

Sun  position  (unit)vector 
Azimuth  angle 

(i  =  1,  2  •  •  •  )  Boundary  values  in  azimuth 
Azimuth  of  target  normal  vector 
Minimum,  maximum  values  of  azimuth  angle 
Nadir  angle 

Nadir  angle  of  target  normal  vector 

Minimum,  maximum  values  of  nadir  angle 

Relative  altitude  parameter,  nadir  angle  of  horizon 

Angular  distance  from  S  to  arbitrary  point  (dAJ  on 
albedo  sphere;  source  to  sun  view  angle 

Angular  distance  from  S  to  R 

Inclination  of  plane  of  Ns  and  R  to  S 

Angle  at  dA2  between  N  and  dAi  ;  target  view  angle 
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<f>e  Angular  distance  of  arbitrary  source  point 

dAj  from  R 

<f>v  Orbital  angular  position,  between  Ns  and  R 

iff  Angle  between  normal  at  d  A4  and  direction  to  dA2 

source  view  angle 
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SECTION  I 
INTRODUCTION 


This  report  extends  one  of  the  problems  discussed  in  an  earlier 
report the  development  of  a  method  for  evaluating  the  "albedo  integral". 
The  aim  of  this  study  is  to  improve  the  speed  at  which  certain  quantities 
are  computed  for  the  control  of  simulation  parameters  in  the  Aerospace 
Environmental  Chamber  (Mark  I).  In  earlier  study  programs,  the  albedo 
integral  was  evaluated  by  strictly  numerical  integration  techniques.  The 
present  seminumerical  method  is  faster,  by  nearly  an  order  of  magnitude, 
than  the  numerical  methods  formerly  used.  This  method  has  been  incor¬ 
porated  into  a  Fortran  language  computer  subroutine. 

A  derivation  of  the  albedo  integral,  for  illumination  intensity,  is 
reproduced  in  Appendix  I,  under  assumptions  that  the  albedo  source  is  a 
homogeneous  sphere  with  a  diffusely  scattering  (Lambert)  surface,  so 
that  the  albedo  is  otherwise  independent  of  surface  and  atmospheric  con¬ 
ditions. 


SECTION  II 

THE  ALBEDO  PROBLEM 


In  order  to  properly  control  the  simulation  of  secondary  radiation 
(albedo  and  planet  radiance)  in  Mark  I,  it  is  necessary  to  determine  the 
illumination  on  an  arbitrarily  oriented  surface  element  at  arbitrary  alti¬ 
tude  and  at  any  position  in  a  trajectory  or  orbital  flight  near  a  reflecting 
celestial  body. 

A  derivation  of  the  albedo  integral,  which  expresses  the  illumina¬ 
tion  intensity,  is  given  in  the  previous  report2  under  assumptions  that 
the  albedo  source  is  a  sphere  having  a  homogeneous,  diffusely  scattering 
surface  so  that  the  albedo  is  otherwise  independent  of  surface  and  atmos¬ 
pheric  conditions.  Then  a  different  primary  body  may  be  distinguished 


■''Cord  H.  Link,  Jr.  "Problems  in  Computing  Radiation  Control 
Functions  for  Mark  I.  "  AEDC-TDR-63 -206,  February  1964. 

2Ibid. 
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by  a  solar  constant  suitable  for  the  distance  from  the  sun,  its  mean 
albedo,  and  its  radius.  This  last  factor  enters  all  secondary  illumina¬ 
tion  calculations  since  they  depend  on  relative  altitude. 


( 0s .  y .  «N  .  /% )  = 


a(<9s  ,y ,  c*n  >^N)  /3(a;0s»y»aN»/^N) 

/  ,/  cos  6  (a,  /3  ;  6S  ,  y)  cos  f  (  a ,  /3  ;  aN  ,  yN )  sin/3d/3d< 

a  (  >  y  >  aN  >  (  a ’ ’  y  ’  aN  >  ) 


The  integral  contains  four  parameters  which  determine  the  con¬ 
figuration  of  boundaries  of  the  surface  over  which  the  integration  is  to 
be  carried  out,  namely,  the  albedo  source  region. 

The  limits  of  integration  are  functions  of  these  same  parameters 
as  well  as  of  the  second  integration  variable.  The  parameters  establish 
the  limits  for  the  second  integration,  as  well  as  controlling  the  func¬ 
tional  form  of  the  integration  limits. 

In  principle  either  of  the  integration  variables  may  be  selected  for 
the  first  integration.  The  azimuth  angle  a  provides  simple  first  integral 
forms,  but  the  function  limits,  involving  the  nadir  angle,  are  sometimes 
double  valued  functions,  a( j8) . 

On  the  other  hand,  the  first  integration  taken  relative  to  the  nadir 
angle  /3  leads  to  more  complex  expressions  for  the  first  integral,  but 
the  functional  limits,  involving  the  azimuth  angle  a,  are  single  valued 
functions,  /3(a). 

This  second  alternative  is  chosen.  Having  once  found  all  the  "anti- 
derivatives"  of  the  integrand  functions  of  j8,  a  simple  differencing  of 
function  values  for  maximum  and  minimum  values  of  /3  (at  a  particular 
value  of  a)  provides  first  definite  integral  numerical  values,  which  are 
now  functions  of  the  four  configuration  parameters  and  a  .  Integration 
over  a  involves  summation  of  first  definite  integral  values. 

Regardless  of  which  variable,  a  or  /3  ,  is  first  used,  when  the  func¬ 
tion  limits  (a({ 3)  or  /3(a))  are  inserted,  some  of  the  expressions 
become  rather  formidable,  and  analytical  evaluations  of  the  second 
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integrals  for  many  of  these  have  not  been  found.  It  is  reasonable  to 
use  simple  numerical  integration  in  the  remaining  variable  a  . 


The  region  of  integration  is  bounded  by  curves  beyond  which  one  or 
more  of  the  integrand  factors  become  negative.  There  are  three  such 
curves  (see  sketch  above).  The  ever-present  horizon  circle  Cy  is  deter¬ 
mined  by  the  relative  altitude  of  the  target  above  the  albedo  source.  The 
terminator  C<9,  the  sunlight- shadow  line,  is  determined  both  by  altitude 
and  by  the  angular  distance  of  the  target  position  from  the  subsolar  point 
on  the  albedo  source.  Finally,  the  ''target  plane  cut"  C£ .  the  intersection 
of  the  plane  of  the  target  with  the  albedo  source,  depends  on  the  specific 
orientation  of  the  target  and  altitude.  The  curve  C$  may  fall  outside 
the  circle,  and  the  curve  C£  does  not  exist  outside  the  horizon  circle. 

So,  depending  on  the  four  parameters,  the  region  of  integration  may  be 
bounded  by  one  curve  Cy,  by  two  curves  (Cy  with  C q,  Co  with  C£,  or 
Cy  with  C£),  or  finally  by  portions  of  all  three  curves. 


Not  only  are  the  boundary  curves  defined  by  the  four  integration 
parameters,  but  their  intersections  are  also,  and  there  may  be  as  many 
as  six  intersections  (see  above).  From  this  arises  part  of  the  complexity 
of  the  problem,  since  the  C£  curve  may  have  any  azimuthal  relation  to  the 
C o  curve,  or  within  the  Cy  circle.  The  logical  sorting  involved  in  deter¬ 
mining  the  boundary  curves  and  their  limits,  for  arbitrary  parameters, 
is  rather  involved  in  the  number  of  decisions  to  be  made.  Yet  for  a 
given  configuration,  only  one  sequence  of  a  few  decisions  serves  to  pro¬ 
vide  all  the  information  required. 
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From  the  standpoint  of  computer  programming,  the  method  described 
here  leads  to  a  large  program,  of  which  only  a  small  part  is  executed  for 
one  given  set  of  parameters.  In  practice  all  the  parameters  may  be  con¬ 
tinually  varying. 

In  the  following  sections,  the  analysis  will  be  developed,  leading  to 
the  computer  program  displayed  herein  as  a  subroutine.  A  logical  flow 
chart  and  Fortran  II  listing  of  the  major  routine  is  given  as  well  as  a 
Fortran  II  listing  of  the  supporting  subroutines.  This  method  turns  out 
to  be  approximately  ten  times  as  fast  in  computing  as  a  corresponding 
purely  numerical  integration  method. 


SECTION  ill 

THE  ALBEDO  INTEGRAL 


The  albedo  integral,  in  its  complete  form,  provides  an  expression 
for  the  intensity  of  illumination  I2  on  an  arbitrarily  oriented  and  positioned 
target  surface  element  A|A2'  attributable  to  albedo  Ae  of  a  homogeneous  dif¬ 
fusely  scattering  sphere  exposed  to  solar  radiation  intensity  Is .  ^  We 
begin  with  the  definitions 


I2 

h  A  e 

/  f  cos  0e  cos  £  sin  /3  d /3 d  a 

(1) 

A  a2 

TT 

cos  0e  = 

COS  ds 

cos  <j6e  +  sin  0S  sin  <f)e  cos  a 

(2) 

cos  £  =  cos 

/3  cos 

/3n  +  sin  /3  sin  /3n  cos  (a  -  ) 

(3) 

sin 

y  =  re/(h  +  re ) 

(4) 

<£e  =  <A  ~  /3 

(5) 

sin  if/  =  sin  /3/sin  y  (6) 


The  integration  is  over  all  a,  /3  within  the  region  where  the  integrand 
factors  are  all  positive.  The  parameters  aN,  /3N  define  the  orientation 
of  a  target  surface  element  to  the  particular  albedo  source  configuration 
defined  by  ds ,  y  • 


3Ibid. 
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For  present  purposes,  the  factor  IsAe  is  taken  as  unity,  leaving  the 
integral 

—  f  f  cos  0e  cos  £  sin  B  d  /3  d  a 

7T 

which  may  be  called  the  "albedo  view  factor".  It  is  a  measure  of  effi¬ 
ciency  of  conversion  of  collimated  illumination  into  scattered  illumina¬ 
tion  on  an  arbitrarily  oriented  surface  element  at  any  point  in  space 
about  a  perfect  diffusely  reflecting  sphere. 

The  integration  relative  to  nadir  angle  /3  is  given  in  Appendix  II. 


a  Geometry  Defining  Location  Parameter  0S(#V/  ^v) 
Fig.  1  Model  Geometry 
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SECTION  IV 

PARAMETERS  OF  ALBEDO  INTEGRAL 


If  S’ is  a  unit  vector  indicating  the  sun  direction,  a  unit  orbital  angular 
momentum  vector,  and  R  the  target  position  vector,  then  a  node  vector 
Ns  may  be  constructed  from  S  x  L  .  Then  the  orbital  angular  position 
4>v  may  be  defined  as  the  angle  between  Ns  and  R  .  The  plane  of  the  orbit 
is  inclined  at  angle  dv  from  the  plane  of  Ns  and  S  .  Then  the  angular 
distance  ds  of  R.  from  S  is  defined  by  (Fig.  la) 

cos  Of,  =  sin  <f>v  cos  dy 
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SECTION  V 

BOUNDARY  CURVES  AND  INTERSECTIONS 


Since  the  boundary  curves  separate  the  region  in  (a,  /3)  for  which  the 
integrand  factors  of  Eq.  (1)  are  positive  from  the  region  where  any 
factor  is  negative,  we  may  write  boundary  equations  as  follows: 

Cy  (Horizon  Circle):  /3  =  y  (7) 

for  all  a . 

C q  (Terminator,  from  Eq.  (2)):  cos  6e  =  0, 
hence 

cos  a  =  -  ctn  0e  ctn  ds  ^g^ 

in  which  Eqs.  (4),  (5),  and  (6)  are  used  to  obtain  expressions  for  a(/3) 
or  /3(a) . 

C£  (Target  Plane  Cut,  from  Eq.  (3)): 

cos  g  =  o, 

hence 


cos  ( a  -  )  =  -  ctn  /3^j  ctn  /3  (9) 

The  intersections  of  Cf  with  Cy  are  obtained  from  Eqs.  (7)  and  (9)  by 
letting  A  a  =  a  -  aN  . 

Then 

cos  A  a  =  —  ctn  /3j^  ctn  y 

N 

ai  =  +  A  a 

k  (MODULO  2n) 

a2  =  aN  -  Act  J  (10) 

The  intersection  of  C q  with  Cy  is  found  by  using  Eqs.  (6),  (5),  and  (7) 
with  Eq.  (8)  as  follows 


P  =  Y 
\}/  -  rr/2 

4>e  =  rt/2  -  y 

(ID 


a ,  =  cos’"1  (—  ctn  ds  tan  y) 
a4  =  2  n  —  a3 
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The  intersections  of  C$  with  C£  are  not  needed  in  the  present  method, 
but  will  be  essential  if  it  is  desired  to  attempt  purely  formal  second-stage 
integration  in  the  future.  The  calculation  of  this  intersection  is  given  in 
Appendix  III;  transformations  between  the  angles  «/r,  t  c£e,  and  B  implied  by 
Eqs.  (5)  and  (6)  are  given  in  Appendix  IV. 


d.  Geometry  Defining  Auxiliary  Angles  (f>e,  ifr 
Fig.  1  Concluded 
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It  is  desirable  to  ensure  that  all  the  angles  of  intersections  (at  ,  a2  , 
a3 ,  and  a*)  of  Eqs.  (10)  and  (11)  are  expressed  as  positive  angles  with¬ 
in  (0,  2 it)  to  eliminate  ambiguities  that  otherwise  occur  in  determining 
when  the  variable  a  is  in  the  range  of  definition  of  one  of  the  curves, 

Co  or  C£. 

It  is  apparent  that  the  limit  points  a,  and  «2  are  symmetrically 
placed  with  respect  to  aN  at  positions  determined  by  y  and  /SN  .  Corre¬ 
spondingly,  a3  and  a4  are  symmetric  relative  to  a  =  0  (located  from 
expressions  involving  y  and  0S).  The  four  integration  parameters  y, 

0S  ,  aN  ,  (SN  remain  arbitrary,  subject  to  limitations 

0  <  y  <  tt/2 
0  <  0S  <77 
0  <  aN  <  2  77 
0  <  /SN  <  77 


SECTION  VI 

MAJOR  DIVISIONS  OF  PARAMETER  RANGES 


All  required  quantities  are  now  defined,  and  we  shall  examine  the 
meaning  of  the  values  of  the  four  parameters  aN  ,  /SN  ,  0S  ,  and  y  .  From 
the  definition  of  y  (Eq.  (4)),  we  find  that  y  approaches  n/2  as  the  altitude 
vanishes,  and  y  approaches  zero  as  altitude  grows  large. 
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We  confine  our  attention  to  the  horizon  circle  and  its  interior, 

/3  <  y  .  In  the  (a,  /3)  coordinate  system  about  the  origin  (the  vehicle  loca¬ 
tion),  a  unit  radius  sphere  is  erected.  The  horizon  circle  is  a  small 
circle,  /3  =  y,  which  entirely  encompasses  the  albedo  source  region,  of 
which  no  part  exists  outside  the  horizon,  /3  >  y  •  The  curve  is  a  great 
circle  on  the  a,  /3  sphere  and  hence  passes  through  the  origin  and  maps 
across  the  interior  of  the  horizon  circle  as  a  straight  line.  The  curve 
C$,  a  great  circle  on  the  albedo  source  sphere,  maps  into  the  horizon 
circle  as  a  part  of  an  ellipse  tangent  to  the  horizon  circle.  The  values 
of  ds  and  /3n  control  the  existence  of  Cfl  and  G£  within  C y  and  the  points 
of  closest  approach  of  these  curves  to  /3  =  o  .  The  missing  ellipse  branch 
of  the  C$  curve,  the  continuation  of  the  terminator,  is  not  defined  in  (a,  /3), 
since  it  is  physically  outside  the  horizon  circle  or  "behind"  it  (see  sketch  * 
on  page  10) , 

Figure  2a  illustrates  the  major  divisions  of  characteristics  imposed 
by  0S  ,  /3jv  ,  and  y  and  by  typical  patterns  of  the  integration  region.  For 
this  illustration,  is  arbitrarily  set  at  77/2  .  Later,  we  shall  examine 
the  influence  of  «N  on  the  problem.  Figure  2a  illustrates  schematically 
some  typical  boundary  patterns. 


0  y  2T  7T  -  y  7r 


a  Typical  Patterns  for  y  =  77/3 

Fig»  2  Major  Divisions  of  Parameter  Ranges  in  the  Horizon  Circle  Cy 
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The  major  ranges  are  noted  in  Fig.  2a  by  use  of  paired  numbers 
(n,,  n2),  the  first  referring  to  the  0S  range,  the  second  to  the  /Sn  range. 

In  range  (0,  0),  only  Cy  bounds  the  region,  and  for  integration  we  have 
0  <  (3  <  y  ,  0  <  a  <  2n  .  This  corresponds  to  a  point  on  the  vehicle 
nearest  the  source  sphere,  near  /3  =  0  and  a  location  of  the  vehicle  not 
far  from  the  subsolar  point  ( 6S  =  0)  on  the  albedo  source  sphere. 

With  no  other  changes,  as  /3n  increases  we  move  from  (0,  0)  to 
(0,1)  where  C£  comes  into  the  horizon  circle.  The  vehicle  itself  begins 
to  mask  part  of  the  source.  The  point  8  =  0  is  still  within  the  source  so 

v 

the  a  limits  are  0  and  277  -,'(3  =  0  is  the  minimum  /3  value  (/3),  and  the 
maximum  (£)  is  either  y  or  dependent  on  a  through  the  equation  for 
curve  Cf . 

When  f3^  =  ff/2.  the  target  plane  (or  its  equivalent  Cf)  bisects  the 
horizon  circle.  Now  a  has  a  range  of  n/2  either  side  of  aN ,  p  is  zero, 
and  /3  is  y  only.  Or  we  may  allow  a  its  full  (0,  277)  range,  but  during 
half  of  this  range  the  curve  Cg  provides  /3  =  0  and  in  the  other  half  Cy 
gives  =  y  while  ^  =  0  . 

As  /3n  grows,  C^moves  into  (0,2),  on  past  the  nadir  point  (3  =  0 , 
and  provides  [3  while  /3  is  y  .  The  range  of  a  is  now  (a,  ,  a 2),  the  region 
where  C£  is  defined.  Finally,  /3 jj  increases  so  far  that  it  is  "on  top" 
of  the  vehicle,  C£  has  swept  completely  across  the  interior  of  the 
horizon  circle,  and  the  integral  value  becomes  zero.  The  target  now 
completely  masks  itself  from  the  albedo  source. 

It  must  be  noted  that  the  integrand  contains  cos  6e,  dc  being  meas¬ 
ured  from  the  subsolar  point  on  the  albedo  sphere.  Thus,  generally, 
the  source  intensity  is  not  symmetric  in  any  way  unless  /3n  =0  or  77  , 
and  these  two  instances  are  not  equivalent  since  one  of  them  includes 
areas  nearer  the  subsolar  point,  and  the  other  is  directly  opposite. 

But  aN  is  arbitrary,  in  practice  a  function  of  /3n  determined  by  the 
vehicle  geometry  and  orientation. 

The  dependence  of  the  Qs  and  (3 N-  ranges  on  the  altitude  parameter,  y, 
is  illustrated  in  Fig.  2b. 

We  now  return  to  case  (0,0)  and  allow  ds  to  vary.  As  we  enter 
(1,0),  the  terminator  Cq  appears  in  the  horizon  circle  at  a  =  77  .  The 
nadir  point  is  still  within  the  integration  region  so  a  ranges  (0,  277), 

/8  =  0  ,  and  /3  are  determined  from  either  Cg  or  C  . 
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Allowing  ds  to  increase  to  n/2  (as  the  vehicle  crosses  over  the 
terminator),  C g  bisects  the  source  field,  and  we  move  on  into  (2,  0). 

V  A 

Here  Co  provides  fi  ,  p  =  y  ,  and  a  ranges  only  over  a3  ,  a4  . 

Finally  0S  increases  so  far  that  C@  leaves  the  horizon  circle  at 
a  =  0,  and  we  have  the  eclipse  condition  where  no  part  of  the  illuminated 
albedo  sphere  is  visible  at  the  vehicle.  The  integral  vanishes.  We  note 
that  cos  £  in  the  integrand  also  destroys  the  apparent  symmetry  in  «  of 
the  source  function  except  in  special  cases.  Instances  where  symmetry 
occurs  were  treated  in  Appendix  II  of  the  previous  report^  as  special 
cases  in  which  the  albedo  integral  can  be  obtained  in  closed  form. 

The  parameter  ranges  labeled  (1,1),  (1,2),  (2,1),  and  (2,  2)  are 
superpositions  of  those  just  described.  The  bounds  are  dependent  on  all 
four  parameters.  When  only  one  of  the  curves  C£  or  C$  establishes  the 
/3 ,  then  the  a  range  necessarily  lies  within  the  corresponding  end  points 
( a,  ,  a2 )  or  ( a3  ,  a4 ) .  More  precise  statements  are  developed  in  subsequent 
sections  as  we  go  more  deeply  into  the  logic  of  sorting  out  the  various 
cases. 


SECTION  VII 

CONFIGURATIONS  OF  ALBEDO  SOURCE  BOUNDARIES 


In  this  section,  we  display  43  distinct  configurations  of  boundaries 
covering  all  useful  values  of  the  four  integration  parameters.  All  of  these 
must  be  examined  for  the  purpose  of  establishing  exact  integration  ranges 
in  a,  ranges  in  which  the  boundaries  are  different  functions  /3(a).  As 
earlier  indicated,  we  shall  eventually  integrate  over  (a,/3)  by  using  exact 
expressions  for  the  first  definite  integral  in  /3  ,  which  contains  functions  /3(a) 
and  /3(a)  ,  then  numerically  integrating  in  a  . 

We  recall  that  all  end  points  (at ,  a2  ,  a3  ,  aA)  are  defined  to  have  values 
in  (0,  2 n)  ,  that  the  curve  C y  is  defined  by  j8  =  y  for  all  a,  that  C£  is 
defined  (Eq.  (8))  in  (at  ,  a2)  ,  and  C$  is  defined  in  (a3  ,  a4)  .  We  do  not  re¬ 
quire  the  intercepts  of  C $  with  C^,  which  would  be  denoted  (a,  ,  a6)  ,  be¬ 
cause  results  based  on  this  knowledge  are  readily  obtainable  by  an  artifice 
which  we  shall  use  in  the  numerical  a  integration.  These  points  would  be 


4Ibid. 
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required  if  one  were  to  attempt  to  find  a  complete  analytical  expression 
for  the  albedo  integral. 


We  return  to  the  notation  of  the  previous  section  for  discussion  of 
the  major  divisions.  Viewed  from  the  origin  of  the  (a,j3)  coordinate 
system,  boundary  curve  Cy  is  a  circle  whose  interior  contains  the 
regions  of  interest.  Curve  C£  is  a  straight  line  segment,  and  Ce  , 
although  an  ellipse  section  tangent  to  C y  ,  is  indicated  as  a  circular  arc 
for  clarity  and  ease  of  drawing.  The  variously  numbered  points  are 
simply  numbered  in  sketches.  The  sides  of  C£  and  Cd  on  which  the  corre¬ 
sponding  integrand  factors  are  positive  are  indicated  by  a  small  arrow, 
pointing  toward  aN  on  C£,  and  toward  a  =  0  on  Co ,  or  to  the  "interior" 
of  these  curves.  Then  the  region  of  interest  is  just  that  part  of  the  pat¬ 
tern  which  is  common  to  the  interiors  of  all  three  curves. 


Case  (0,  0)=  pN  <2-y,  0$  <y 
0  <  a  <  27r ,  (3  3  0  ,  (3hy 


Case  (1, 0):  pN<|-y,  y  <  9S<| 

0  <  a  <  27r  p  =  0  p  =  y 

or  p  *  p  (Cq)  for  aj  <  a  < 


Case  (2,0):  PN  <  |-y  l  <  0$  <  tt  -  y 


&  <  a  <  CI3 
and  <  a  <  2 7r 


split 

scan 


p  =  p(C9)  p  =  y 


Case  (0, 1):  9S  <  y 
0  <  a  <  2  7T 

v  "  ~ 

p  =  0 

/\ 

p  =  y 

or  p  =  p(C^) 

if  I  °<a  <a2  | 

/  a  i  <  a  <  2  7r ) 


split 

test 


0  <  a  <  2  7r 

r-o" 

B  •  y 
or 

P  ■  P<C^> 
if 

aj  <  a  <  d2 


N 


Case  (0, 2):  0S  <  y 

P  =  y 

B  ■  P<Cg) 
a2  1  a  1  al 


f<PN<f +  r 


P  =  y 

B  =  P  <Cg) 

0  <  a  <  oq  | 
and  a2  <  a  <  2?r  I 


split 

scan 


Fig.  3  Configurations  of  One  and  Two  Boundary  Curves 
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Cases  (o,0),  (0,1),  (0,2),  (1,  0),  and  (2,  0),  shown  in  Fig.  3,  are 
largely  self-explanatory.  Cases  (0,1),  (0,2),  and  (2,0),  however,  intro¬ 
duce  the  problem  of  the  "split  range".  Although  the  two  orientations  shown 
in  (0,1)  and  (0,2)  are  geometrically  equivalent,  they  are  logically  distinct 
since  all  intersection  values  are  defined  in  (0,2n-)  .  For  example,  in  case 
(1,0)  when  the  order  of  end  points  is  a,  <  a2  ,  the  curve  is  defined  in 
(a,  ,  Oj)  ,  but  when  the  order  is  a2  <  a,  ,  is  defined  in  the  split  range 
(0,  aj)  and  (a,  ,  2n)  .  Thus,  the  order  of  end  points  is  essential  to  the 
orderly  determination  of  the  range  in  which  a  may  be  during  numerical 
integration  and  for  the  selection  of  the  proper  boundaries  (^3  ,  /3)  for  a 
given  a . 

The  split  range  is  also  used  in  the  initiation  and  advance  of  a  during 
numerical  integration.  If  the  range  is  split,  a  is  scanned  over  the  two 
parts  successively. 

Case  (1,  1)  is  the  most  complicated  group  of  configurations  because 
the  a  range  is  0,  2n  and  both  C£  and  C$  are  present.  Fig.  4.  Each  con¬ 
figuration  is  labeled  by  letter  referring  to  the  corresponding  permutation 


Three  Boundary  Curves 

f-r<  PN<|,  r<es<?,  P  -o  ,  0<a< 2u 


Configuration 

a  Order 

0oundary  Curve  for  p 

Notes 

A 

(0)  1  2  3  4  (2fl)  y  4 

r  0 

r 

0 

1324 

r  4 

(0, 4)  0 

r 

C 

1342 

Y  4 

4  4 

r 

D 

2314 

4  r 

0  (0,4) 

4 

E 

23  4  1 

4  r 

0  y 

4  - 

F 

3124 

r  0 

(0,4)  0 

r 

G 

3  142 

r  0 

(0,4)  4 

r 

H 

3  24  1 

4  (0,4) 

0  y 

4  * 

1 

3412 

r  0 

r  4 

r 

Symbol 

Meaning 

(0,4) 

p  =  Min(p(C0),  p(C£)) 

r 

P  =  r 

e 

t  3  P(c0) 

4 

p  -  p(C|) 

* 

Split  Test 

Fig, 

,  4  Case 

(i,  i) 
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of  end  points,  given  in  the  table  at  the  bottom  of  the  figure.  Beside  each 
permutation  appears  the  order  of  subscripts  in  the  boundary  curves  from 
which  fi  is  to  be  found.  For  this  case,  the  point  /S  =  0 ,  lying  interior  to 
the  integration  region,  is  also  /3  =  0  .  Where  a  double  subscript  occurs, 
we  make  use  of  the  artifice  (previously  mentioned)  to  determine  whether 
to  use  Cq  or  C^-.  Here,  when  /§  is  to  be  found  in  a  range  of  a  where  both 
C,9  and  are  defined,  we  use  the  C#  and  Cg  definitions  to  determine 
both  values  of  /3  ,  i.  e.  ,  B{Co)  and  /3(Cfi)  ;  we  then  select  the  least  of  these 
to  be  /S  . 

We  give  an  example  of  interpretation  of  the  table  of  Fig.  4.  Select 
configuration  B.  Initiate  a  at  0  ,  change  the  value  by  the  (fixed)  step 
size,  and  test  a  to  see  when  it  is  in  each  subsequent  range.  In  (0,  a3)  , 

@  =  y  ;  in  (a,  a3 )  ,  /3  =  /3  (C£)  .  In  ( a3  ,  ctj )  ,  /3  min  [  /3  (C$) ,  /3(Cg)  \  ;  then 
in  ( a2  ,  a„),  jS  =  (3{Co)  .  Finally,  in  («4  ,  277),  /3  =  y  . 

Note  that  in  any  configuration  in  which  a2  <  a,  ,  the  range  test  must  be 
split,  as  noted  earlier. 

V 

Case  (1,  2)  in  Fig.  5  has  /3  defined  by  the  curve  C^,  and,  hence,  a  is 
scanned  only  over  the  range  of  definition  of  C^;  that  is,  over  (a3  ,  a2)  if 
a,  <  ctj  ,  or  over  (0,  a,)  and  (a,  ,  2tt)  if  a,  <  a,  .  In  the  table,  the  symbol  [&] 

A  V 

means  that  we  use  /3  =  /3(C<?)  only  if  it  is  greater  than  /3  =  /9(C£)  ;  otherwise, 
there  is  no  contribution  to  the  integral  for  the  current  value  of  a. 

V 

Case  (2,  1)  in  Fig.  6  has  /3  defined  by  C$ ,  so  a  is  scanned  over  (0,  a3) 
and  ( at ,  2n) ,  a  split  scan.  The  notation  <f*  in  the  table  means  that 
=  /3(C£)  if  only  /§  >  ft  ;  otherwise,  the  current  value  at  a  contributes 
nothing  to  the  integral. 

A  V 

Finally,  case  (2,  2)  in  Fig.  7  has  (3  defined  by  =  y  ,  and  we  select  '/3 
as  the  greatest  of  /3( C^)  and  /3(C#)  ,  which  is  the  meaning  of  the  symbol 
(0,  £)  in  the  table.  Note  that  configuration  G  illustrates  that  nonover¬ 
lapping  of  the  ranges  of  (a3  ,  a2)  and  («3  ,  a4)  leads  to  zero  value  of  the 
integral.  For  case  (2,  2)  we  let  a  scan  only  the  least  of  the  spans  of  C 0  or 
C£ ;  if  this  is  C<?  then  the  a  scan  is  split,  but  if  C£  ,  the  a  scan  may  or  may 
not  be  split.  Notations  of  split  scan  test  appear  on  the  figures. 

This  completes  the  details  of  the  logical  procedures  for  doing  the 
numerical  integration  in  a  .  From  the  tables  on  the  figures,  the  logic  flow 
chart  in  Appendix  V  was  derived;  the  problem  was  then  programmed  for 
computer  directly  from  the  flow  chart.  The  Fortran  listing  is  shown  in 
Appendix  VI. 
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Three  Boundary  Curves 

f<pN<f  +  r,  r<es<| ,  K-PKj) 

if  a2  <alf  then  a2<  cKc^  ;  but  if  ax  <  a2  , 
then  (KcKaj  and  a2  <  a  S  2?r 


Configuration  a  Order  Boundary  Curves  for  p  Notes 

A  (0)  1  3  2  4  (2tt) 

B  1342 

C  2134 

D  23  14 

E  23  4  1 

F  3142 

G  3214 

H  3  24  1 

|  3  4  2  1 


Three  Boundary  Curves 

f-r<PN<f  .  f<9s<rr-  y  ,  p  =  p(C9) 


Split  Scan  0<a<a3  and  a4<a<27r 


ABC 


Configuration 

a  Order 

A 

Boundary  Curves  for  p 

Notes 

A 

(0)  1  23  4  (271) 

r  |*  r 

- 

y 

* 

B 

1324 

y  - 

- 

y 

0 

C 

2314 

|*  r  - 

- 

t* 

D 

23  4  1 

i*  y  - 

y 

* 

E 

3124 

r  -  - 

- 

y 

* 

F 

3  142 

r  -  - 

i* 

y 

■> 

G 

3214 

-  - 

- 

- 

** 

H 

3  24  1 

£♦  -  - 

y 

i* 

1 

3412 

y  -  y 

V 

y 

0 

y  -  -  00]  y 

r  -  -  -  r 

-  y  -  -  - 

-  y  ISO  -  - 

-  y  0  y  - 

y  [0]  -  -  y 

-  -  09]  -  - 

-  -  [B1  y  - 

-  -  -  r  - 


Symbol  Meaning 

[0]  $  -  pfCg)  only  if  P  >  P,  Otherwise 

No  Contribution  to  integral 

Means  No  Contribution  to  Integral 
*  Split  Scan  in  a 


Symbol  Meaning 

C  P  -  P(C^)  if  p  >  P 

Otherwise  No  Contribution 

*  Split  Scan 

*«■  Identically  Zero  Integral 


Fig,  5  Case  (1,  2) 


Fig.  6  Case  (2,1) 
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Three  Boundary  Curves  +  y  ,  |<0s<7r-y,p=y 


a  Ranges  over  Least  Span,  Covering  Range  of  C|  if  A  a  <  a 3, 


Otherwise  over  Cg.  Span  of  Cg  Always  Split  0  <  a  <  03  and 


Configuration 

a  Order  Boundary  Curves  for  P 

Notes 

A 

(0)  1  3  2  4  (2tt)  (9,1)  -  -  - 

(9, 1) 

* 

B 

1342 

(0, 4)  -  -  - 

(0, 1) 

C 

2134 

-  (0  1)  -  - 

- 

D 

23  14 

-  (9,  £ )  -  - 

- 

E 

3  124 

0  -  - 

0 

* 

F 

3  142 

(0,|)  -  -  - 

<e,  |) 

* 

G 

3214 

- 

- 

H 

3  2  4  1 

-  -  (0,|) 

- 

1 

3  4  2  1 

-  -  (0,|) 

- 

Symbol 

Meaning 

(0,  £ ) 

p  =  Max  [P(Cg),  p(C|3 

❖ 

Split  Test  if  a  Ranges  over  C| 

** 

Identically  Zero  Integral 

Fig.  7  Case  (2,  2) 


19 


A  E  DC-T  R-65-202 


SECTION  VIII 
CONCLUSIONS 


In  an  attempt  to  gain  computing  speed  in  the  evaluation  of  the  albedo 
integral,  a  double  integral,  the  problem  has  been  changed  from  a  straight¬ 
forward  numerical  integration  to  a  much  faster  but  more  complex  semi- 
numerical  integration.  For  example,  whereas  formerly  the  a,  fi  range 
was  covered  by  a  mesh  of  72  X  36  points,  the  present  method  requires 
somewhat  more  computation  per  point  through  a  more  complex  logic  net¬ 
work  at  only  72  points,  and  the  accuracy  is  improved  by  the  formal  first 
integration.  The  time  improvement  is  approximately  one  order  of 
magnitude. 

In  application,  a  particular  Mark  I  control  program  formerly  required 
from  25  to  35  min  (IBM  7074)  to  generate  simulation  parameters  for  a 
90 -min  orbit  with  a  simulation  interval  of  two  minutes.  The  same  results 
are  now  produced  in  approximately  four  minutes. 

It  appears  unlikely  that  significant  gains  in  computing  speeds  can  be 
made  by  using  a  purely  formal  solution  to  this  problem.  Second  integrals 
will  contain  many  more  terms,  some  quite  complex,  and  much  of  the  gains 
made  by  having  a  single  evaluation  to  perform  will  be  lost  in  the  sheer 
bulk  of  the  expressions  involved.  Most  of  the  logic  of  the  present  method 
would  still  apply  for  selecting  integration  limits  and  function  groups  to  be 
evaluated.  Some  gain  may  result  in  changing  the  variable  of  first  integra¬ 
tion,  and  this  will  be  studied  in  the  future.  It  may  also  be  possible  to 
develop  rapidly  computing  empirical  approximating  functions,  especially 
over  limited  ranges  of  the  integral  parameters. 
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APPENDIX  I 

DERIVATION  OF  THE  ALBEDO  INTEGRAL 


In  the  following  discussion,  the  albedo  source  body  is  taken  to  be 
the  earth.  Substitution  of  appropriate  values  for  radius,  albedo,  and 
solar  constant  allows  extension  to  any  source  body. 

To  compute  earth  albedo  and  radiance  integrals  for  a  surface  element 
having  arbitrary  orientation  and  position,  we  make  the  following  assump¬ 
tions: 

1.  that  Albedo  is  a  uniform  property  of  the  earth's  surface, 

2.  that  the  earth  is  a  sphere, 

3.  that  the  earth's  surface  is  diffusely  reflecting,  and 

4.  that  the  earth  has  no  atmosphere. 

Is  Solar  constant,  intensity  of  solar  radiation  at  earth 

Ae  Albedo,  fraction  of  solar  constant  reflected,  a 
surface  property. 

The  solar  radiation  incident  on  an  area  element  AAX  having  its  normal 
inclined  at  angle  0e  to  sun  direction  is 

f  ls  COS  0e  AAX  f°r  cos  -  0 

AIe  =  1 

l  0  for  cos  9e  <  0 

Of  this  a  fraction  Ae  is  reflected  diffusely  by  AAX ;  hence,  the  intensity 
per  unit  solid  angle  AI^  in  a  direction  inclined  at  angle  if/  to  the  surface 
normal  is 

Ae  Is 

AI ^  =  — - cos  6e  cos  if/  AAX 

The  intensity  included  in  solid  angle  A<d  is 

-A  I 

A I6J  =  — L — L_  cos  0e  qos  if/  A  Ax  A<u 

An  area  element  AA2,  at  distance  pe  from  AAX,  having  its  normal 
inclined  at  angle  £  to  the  direction  of  pe  ,  intercepts  a  solid  angle 
(Fig.  1-1). 

.  A  A,  cos  A 

A  co  = - — = - - - 

Pe 
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Hence  the  intensity  arriving  at  AA2  is 


AI2 


K 

7T 


I. 


cos  de  cos  \p  cos  f  AAt  A  Aj 


Let  the  area  element  AA2  be  located  at  altitude  h  above  earth  of 
radius  re  .  From  this  point,  the  portion  of  the  earth  that  can  be  seen  is 
confined  within  a  horizon  circle.  The  angle  y  between  the  direction  of 
earth  center  and  the  horizon  circle  is  defined  by 

sin  y  =  re/(re  +  h )  (0<y<  — ) 

At  the  earth  center,  let  6S  be  the  angle  between  the  direction  to 
AA2  and  the  sun  direction;  let'  de  be  the  angle  between  the  area  element 
AA2  on  earth  and  the  sun;  let  <£e  be  the  angle  between  the  directions  of 
AAj  and  AA2 . 
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x  =  re  sin  i|i  =  (re  +  h)  sin  p 
Horizon  Line  s'n  ^  re/(re  +  h) 


T  sin  y 
i(i  =  3  +  0 


To  Sun 


Earth  Surface 
Radius  rc 


sin  y  ■  re/(re  +  h) 

cos  8e  =cos  0  cos  0S  +  sin  0  sin  6S  cos  a 

sin  =  sin  p/sin  y 

V  "P  +  0e 

Fig.  1-2  Coordinate  Relations  for  Albedo  and  Earth  Radiance  Calculation 


At  the  element  AA2  .  if/  is  the  angle  between  the  normal  to  At  and  the 
direction  of  AA2  as  before. 

At  the  element  AA2 ,  let  fi  be  the  angle  between  the  direction  to  earth 
center  and  AA2 .  Angle  /3  is  a  "nadir"  angle. 

Then  the  following  relations  hold  (Fig.  1-2): 

'A  =  Z3  +  4> 

sin  ip  =  sin  /3  /  sin  y 

cos  6e  =  cos  9S  cos  <A  +  sin  9S  sin0e  cos  a 

where  a  is  the  angle  about  the  line  from  AA2  to  earth  center  measured 
from  the  plane  including  this  line  and  the  sun.  Angle  a  is  the  azimuth 
angle  (Fig.  1-3). 
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Now  the  element  AA,  is  described  in  spherical  coordinates  having 
polar  axis  along  the  earth-to-  AA2  line,  longitude  a,  and  co-latitude  <fie: 

AAi  =  re2  sin  <pe  d  a  d  <pe 


The  a,  3  Coordinate  System  at  AA2 


P  -0 

cos  £,  =cosPco$PN  +  sinPsinPNcos(a-aN) 


Fig.  1-3  Albedo-Radiance  Integration  Coordinates 


Similarly,  we  may  describe  a  spherical  area  element  in  terms  of 
a ,  /3  ,  and  pe  from  AA2  : 

AA  p  =  p*  sin  £  d/3  da 


Hence  any  point  may  be  described  by  either  (re,  cpe,  a)  or  (pe  ,  (3  ,  a) 
at  AAt  ,  and  we  find  that  AAp  and  AA!  are  related  by  simple  projective 
properties: 


so  that 


AAt  cos  ^ 


A  A  =  Pe2  sin  fi  d  /3  d  a 

1  COS  Ip 
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We  may  now  write  the  intensity  at  AAa  caused  by  reflection  from  AAt  : 


A  l_2 
A  A, 


A  e  I. 


A  I 

e  1  s 


cos  0C  cos  \jj 


('os  f  pe2  sin  j 


ios  \Jj 


d/3 


:os  6e  c £»s  f  sin  P  d/3  d  a 


Integrating  over  the  interior  of  the  horizon  circle,  we  obtain 

1  AI  C  '  »amax 

-,-2 —  ( A  Aa ,  ,  y ,  £ )  =  — 2 — —  J  /  cos  6e  (d. ,  j8  ,  0S ,  y )  cos  £  sin  jS  d  j8  d  a 

AAa  -  77  B  =  0  Ja  ■ 

r*  v  in 

where  am[n  ,  amax  may  be  functions  of  /3 ,  and  there  may  be  several  distinct 
regions  or  functions. 


In  the  spherical  coordinates  a ,  /3  ,  any  orientation  of  surface  ele¬ 
ment  may  be  described  by  the  direction  angles  of  its  normal,  /3N. 

Then  the  angle  £  is  found  from 

cos  £  ~  cos  jS  cos  /3n  4-  sin  /3  sin  /3ft  cos  (a  —  ) 

Since  0e  =  \fj  —  (3 ,  and  sin  =  sin  /3/ sin  y 

cos  0e  =  sin  ip  (cos  0S  sin  /3  4-  sin  0S  cos  /3  cos  a) 

4-  cos  ijj  (cos  6S  cos  /3  —  sin  0S  sin  /3  cos  a) 


The  intensity  integrand  is  completely  expressible  in  the  two  variables 
a  ,  /3  ,  and  the  configuration  parameters  0S,  y  ,  ajj  ,  /3 n  • 

At  this  stage  it  is  possible  to  integrate  numerically  by  letting  a  range 
from  0  to  2 77  and  /3  range  from  0  to  y,  provided  that 

cos  £  >  o 
cos  6>e  >  0 
cos  ifj  >  0 

and  using  (=0)  for  any  («,/3)  violating  these  conditions. 

Ignoring  for  the  moment  the  constant  Aie‘Is/;r,  we  have  the  following 
terms  to  be  integrated  over  a  and  /3  : 


1. 

cos 

Pn 

cos 

% 

cos2  /3 

sin 

P 

cos 

0 

2. 

—  cos 

/3n 

sin 

6 

s 

cos  /3 

sin2 

P 

cos 

ifj  cos  a 

3. 

cos 

/3n 

cos 

6 

s 

cos  /3 

sin2 

P 

sin 

* 

4. 

cos 

sin 

(9 

s 

COS2  /3 

sin 

P 

sin 

ijj  cos  a 
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5. 

sin 

/3n 

cos 

e 

s 

cos 

«N 

sin2 

ft 

cos 

ft 

cos 

4> 

COS 

a 

6. 

sin 

/3  N 

cos 

e 

s 

sin 

«N 

sin2 

ft 

cos 

ft 

cos 

sin 

a 

7. 

—  sin 

/3n 

sin 

9 

s 

cos 

«N 

sin3 

ft 

cos 

cos2 

a 

8. 

-  sin 

/3n 

sin 

9 

s 

sin 

«N 

sin3 

ft 

cos 

sin 

a 

9. 

sin 

/3n 

cos 

9 

s 

cos 

«N 

sin3 

ft 

sin 

COS 

a 

10. 

sin 

/3n 

cos 

9 

s 

sin 

«N 

sin3 

ft 

sin 

sin 

a 

11. 

sin 

sin 

9 

s 

cos 

«N 

sin2 

ft 

cos 

ft 

sin 

cos2 

a 

12. 

sin 

/3n 

sin 

9 

s 

sin 

«N 

sin2 

ft 

cos 

ft 

sin 

* 

sin 

a 

Note:  sin  i/j  =  sin  /3/sin  y 

As  long  as  there  are  no  boundaries  of  earth  surface  for  which 
a  =  a  ( /3  ),so  that  a  can  range  from  0  to  2  77 ,  we  may  integrate  relative  to 
a  and  obtain  simple  results.  Integrals  (1,  3)  do  not  contain  a  so  the 
integration  results  in  a  factor  2n .  Integrals  (2,  4,  5,  6,  8,  9,  10,  12) 
contain  only  sin  a  or  cos  a  and  vanish.  Integrals  (7,  11)  contain  cos2  a 
or  cos  a  sin  a  and  result  in  a  factor  of  tj  . 


These  conditions  are  satisfied  as  long  as  we  have  both 

ft  N  -  ~y~  -  y 

0S  <  y 

If  >  tt/2  +  y  or  9s  >  n  -  y,  the  earlier  conditions  on  cos  £  or 
cos  9e  are  violated  and  the  entire  integral  (I2/AA2)  vanishes. 

For  -j-  -  y  <  <  -2-  +  y 

and/  or  y  <  0g  <  77  -  y 

there  exist  boundaries  of  form  a(/3),  and  the  integration  becomes  com¬ 
plicated.  The  integration  is  bounded  by  arcs  of  one,  two,  or  three 
curves  of  a(/3),  whose  intersections  are  generally  given  by  implicit 
functions.  A  first  integration  may  be  done  formally;  expressions  result 
for  which  the  integrals  are  not  available  in  closed  form. 


Numerical  integration  may  be  accomplished  in  an  easily  compre¬ 
hended  manner  by  referring  to  the  earlier  integral  expression.  The 
product  cos  f  cos  0p  sin  j3  may  be  calculated  term  by  term,  and  in  addi¬ 
tion,  the  expression  cos  ifr  can  be  evaluated  to  ensure  that  the  conditions 


cos  £ 

cos  0 
e 


>  0 


cos  \fj  J 

are  satisfied.  For  some  value  combinations  of 
tions  can  be  carried  out. 


aN,  /3„,  <9b  the  integra- 
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APPENDIX  II 

FIRST  INTEGRATION  IN  NADIR  ANGLE  /3 


Table  II -I  displays  the  twelve  possible  integrands,  with  their  param¬ 
eter  coefficients.  The  last  eight  of  these  may  be  grouped  in  pairs  and 
combined  by  use  of  the  identity 

cos  a  cos  a ^  +  sin  a  sin  =  cos  (a  —  aj^) 

Further  grouping  can  then  be  performed  based  on  the  formal  similarity 
of  integrands.  We  use  the  numbering  of  Table  II -I  to  identify  integrands 
and  the  following  definitions  of  parameter  functions. 


Ai 

=  sin  /3n  sin  (9g  cos  a  cos  (a  - 

aN ) 

a2 

=  COS  /3^  COS  0S 

a3 

=  sin  cos  0Q  cos  ( a  -  ) 

Ah 

=  cos  j3^  sin  (9s  cos  a 

(1) 

a2 

j *  cos2 /3  sin  /3  (1  —  sin2  /3  /  sin 

;  d/3 

(2) 

-A4 

J *  cos  /3  sin2  jS  (1  -  sin2 /3 /sin2  y 

;  d/1 

(3) 

a2 

j *  cos/3  (sin3  /3/sin  y)  d/3 

(4) 

Ah 

J *  cos2  /3  (sin2  /3/sin  y )  d /3 

(5,  6) 

a3 

J  cos  /3  sin2  /3  (1  —  sin2/3/sin 

d/3 

(7,  8) 

-A, 

^  sin3  ^  (1  -  sin2 /3 /sin2  y)^ 

d/3 

(9,  10) 

A3 

j '  (sin4  ^8/ sin  y)  d /3 

(11,  12) 

a. 

f  cos  /3  (sin3  /3  /.sin  y )  d /3 
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TABLE  ll-l 

INTEGRAND  FORMS  FOR  ALBEDO 

1.  cos  P|\j  cos  0S  cos2  p  sin  p  cos  V 

2.  -cos  pN  sin  0$  cos  p  sin2  p  cos  cos  a 

3.  cos  pj^  cos  0$  cos  p  sin2p  sin  ip 

4.  cos  PN  sin  0$  cos2p  sin  p  sin  cos  a 

5.  sin  pN  cos  0S  cos  aN  sin2  p  cos  p  cos  cos  a 

6.  sin  PN  cos  0$  sin  aN  sin2  p  cos  p  cos  V  sin  a 

7.  -sin  PN  sin  0$  cos  aN  sin^  p  cosy  cos2 a 

8.  -sin  P|\|  sin  0$  sin  sin^  p  cos  y  sin  a  cos  a 

9.  sin  PN  cos  0$  cos  aN  sin^  p  sin  ip  cos  a 

10.  sin  pN  cos  0  sin  aN  sin^  p  sin  y  sin  a 

11.  sin  PN  sin  0$  cos  aN  sin  p  cos  p  sin  y  cosc  a 

12.  sin  P|\j  sin  0$  sin  sin2  p  cos  P  sin  vp  sin  a  cos  a 

Note:  sin  y  =  sin  p  / sin  y 
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Regrouping  for  formal  similarity 

(1,  7,  8)  (A2  +  A,)  f  cos2  /3  (1  -  s  in2  /3  /sin2  yY2  sin  /3  d/3 

-  Aj  J  (1  -  sin2  /3/sin2  y  )  ^  sin /3  d/3 

(2,5,6)  ( A3  —  A4 )  J"  sin2 /3  (l  —  sin2  /3 /sin2  y)  ^  cos /3  d /3 

(3,  11,  12)  (At  +  A2 )  J  (sin3 /3/sin  y)  cos/3  d/3 

(4,  9,  10)  A4  £  (sin2  /3 /sin  y )  d /3  +  ( A3  -  A4 )  J"  (sin*  /3 / sin  y)  d /3 

We  proceed  to  integrate,  first  making  the  following  substitutions 
(1,  7,  8)  Let 


sin  y  =  G 

cos  y  *=  B  cos  /3  ~  x 

and  note  that 

.  2 

sin  y  - 

sin2  /3  *»  cos2  /3  —  cos2  y 

sin  (3  d (3  =  -  d  cos  (3  =  ~  d  x 

Then  we  obtain 

(  Aj  +  Aj  ) 

G 

/*  x2(x2  -  B2)I/2 

dx  +  f  U  -  B*Y2  dx 

(  A2  +  Ax ) 

G 

[f  <*'  - 

‘  +  —  X  (x2  -  B2  Y2  ~  ~  In  (x  +  (x2 
8  8 

-  B2  )Va )  ] 

f-i  (x2  -  B2),/2  -  Zn  (x  +  (x2 

-  B2)Va)] 

( A2  +  Al ) 

G 

4-  (x2  -  b2)3/2 

4 

1 

f  Ai  Ax  +  A2 

— |x  (x2  -  B2)^  -  B 2  In  (x  +  (x2  - 

-  B2 ) 1/2 )  } 

T 

L  2G  G' 

Finally 

(1,  7,  8)  ^-{[2A,  -  (Aj  +  A2)  B2]  [x  (x2  -  B2)^  -  B2  In  (x  +  (x2  -  B2 )^  ) J 

-  (A,  +  A2)  x  (x2  -  B2 ) 34 } 
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(2,  5,  6)  Let 


where 


sin  y  —  G  sin  /3  =  y  cos  if/  =  z 


sin  xjj  =  sin  /3 / sin  y 


and  note 


cos  /3  d/3  =  d  sin  /3  =  dy 


Then  we  obtain 

(  A3  2  /p2  2,/2  J 

- £ -  J  y  (G  -  y  )  2  dy 


( A3  —  A, 


G 

a4  -a3 

4  G 

a4~a3 

4  G 


~  [~L(°2  “  y2)i/2  +  "1"  {y(G’“  y2 A2  +  G2  sin_1 

[y(G2  -  /)Vs  -  L  | y  ( G2  -  y2)1/2  +  G2  t/rjj 


-i  y_ 
G 


yz3-  -i-(yz+Gi/f) 


(3,  11,  12)  Let 

then 

(A, 


(4,  9,  10)  Let 
then 


G  =  sin  y 


sin3  /3  cos  ft  d  ft 
sin  y 


(  A2  +  Aj  ) 
4  G 


sin*4  /3 


G  =  sin  y 


f  sin 2  ft  1ft  +  (^3~g-4-  /  sin4  ft  d/3 

=  Af  [y  (£  “  sin^  cos  /3)j 

( A 3  “  A4 )  f  sin3  B  cos  B  3  [1  0  ^  vl 

+  G  |_ - ^ +  T  IT  “  sin/3  008 

=  — 7T  f  (A<  -  A3 )  sin3  (6  cos  /3 


+  -j-  (A*  +  3A3)  (ft  -  sin  ft  cos  ft) 
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Now  all  groups  have  a  common  factor  1/(4G);  this  factor  is  ignored 
in  practice  until  final  calculation  of  the  albedo  view  factor,  which  is 
the  calculated  value  of  the  integral  multiplied  by 

1/(4  n  sin  y) 

The  actual  albedo  illumination  intensity  is  then  gotten  by  multiplying 
by  the  albedo  Ae  and  solar  constant  Is  • 
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APPENDIX  III 

INTERIOR  INTERSECTIONS  OF  BOUNDARY  CURVES 


Intercepts  of  C$  with  C£  are  defined  by  the  system  of  equations 

cos  a  =  —  l/(tan  6S  tan  (j>e )  (III-l) 

(III  -  2 ) 


with  the  conditions 


cos(a-aN)  =  -  l/(tan  /3N  tan/3) 


4>e  =  'A  “  P 

sin  iff  =  sin  /3 /sin  y 


(III  — 3 ) 
(III-4) 


thus  all  four  parameters  6S  ,  y,  aN ,  and  are  involved.  From  the  condi¬ 
tions  of  Eqs.  (Ill- 3 )  and  (III-4)  we  derive  the  relation 


tan  /3 


0 


sin  cb  sm 


y 


1  —  sin  y  cos  <f> 


(III- 5) 


We  expand  the  left  side  of  Eq.  (Ill- 2 )  and  use  Eq.  (Ill- 5)  on  the  right  of 
Eq.  (III-2)  to  obtain 


where 


Let 


C  cos  a  +  S  sin  a  =  (G  cos  (f>e  -  1)/BG  sin  (f>e 
B  =  tan  /SN 

C  =  cos  ,  G  =  sin  y  ,  S  =  sin 
T  =  tan  6S 


(III- 6 ) 


and  substitute  Eq.  (Ill- 1)  on  the  left  of  Eq.  (Ill- 6)  to  obtain 


+  S  1  - 


C  cos 
T  sin 

Now  write  the  middle  term  as 

s 


3°S  ( A 


T 2  sin 2  <f> 


% 


T  sin  0e 


(T2  sin2  <f>  -  cos2  </>  )  ^ 


T  sin  <fi 


G  c  os  —  1 
B  G  sin  cf>e 


T2  -  (T2  +  1) 


4>e] 


and  factor  out  sin  <f>e  assuming  <£e  ^  0.  Isolate  the  radical  on  the  left  side 
and  obtain,  by  squaring  and  rearranging, 

G2  [(B2  +1)(T2  +1)  -  (1  -BCT)2]  cos  <£e  -  2T  (BC  +T)  G  cos  <f>e  +  T2  (l  -S2B2G2)  =  0 

(HI- 7) 
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Solving  this  quadratic  for  cos  <f>e,  we  find 

.  ,  =  _T_  _  (BC  +  T)  ±  {  (BC  +  T)2  ~  [  (B2  +  1)  (T2  +  l)  -  (  1  -  BCT  )2  ]  (  1  -  S2  B2  G2 )  )  \ 

6  G  (  B2  +  1  )  (T2  +  1 )  -  (  1  -  BCT)2 

(HI- 8) 

From  these  two  (±)  values  of  cos  <f>e  we  may  find  corresponding  values  for 
sin  (f>e,  noting  that  <£e  <  -| —  y  always,  by  definition.  Use  the  values  of 
cos  <f>e,  \  sin  <f)e  in  Eq.  (Ill - 5)  to  obtain  two  corresponding  values  of  tan  f3  , 
hence  fi  ,  noting  that  for  this  purpose  (3  <  y  .  At  the  same  time,  Eq.  (Ill - 1) 
allows  evaluation  of  the  two  values  of  a  ,  and  the  intersection  points  for 
C o  and  can  be  found.  These  points  are  labeled  (a5  ,  /3S )  and  (a6  ,  /S6)  . 

Now  the  discriminant  of  Eq.  (Ill- 8)  provides  indications  of  the 
presence  of  two,  one,  or  no  intersections  of  Co  and  . 

From  the  fact  that  (a5  ,  a6)  contains  a  values  of  intersection  for  two 
curves  defined  between  both  pairs  (a,  ,  a2)  and  (a3  ,  a4)  ,  of  necessity 

( as  ,  a6 )  values  lie  inside  both  ranges  (a,  ,  Oj)  and  (a3  ,  aA) .  That  is, 

( as  ,  a6 )  values  are  within  the  a  ranges  of  definition  which  are  common  to 
both  Cf  and  Cq  •  This  is  the  region  where,  in  practice,  we  may  compute 
both  /3( Cq)  and  /3(C|)  to  select  which  shall  be  used.  In  this  way  we  avoid 
calculation  of  ( a5  ,  a6)  and  avoid  further  complication  of  the  selection 
logic. 
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APPENDIX  IV 

TRANSFORMATIONS  FOR/3,0e,  AND  0 


Transformations  for  xp,  0e,  and  /3  are  based  on  the  relations 

0  =  0e  +  Z3 

sin  ip  =  sin  /3  / sin  y 

.2  .  V 

F  =  (1  +  sin  y  —  2  sin  y  cos  0e)  2 

The  transformations  are 

sin  xp  =  sin /3/sin  y  =  sin0e/F 

cos  xp  =  (sin2  y  -  sin2  j3)^/ sin  y 


=  (cos  0e  -  sin  y)/F 


sin  j8  =  sin  y  sin  0/F  -  sin  y  sin  0 
cos  j8 
sin  0e 


(1  —  sin  y  cos  0)/F  =  (1  —  sin2  y  sin2  0)^ 


|^cos  j8  -  (sin2  y  -  sin2  )  2  j  sin  0  /  sin  y 
=  sin  \p  £(1  —  sin2  y  sin2  xp)  2  —  sin  y  cos  xp  J 

cos  <pe  =  ^sin2  /3  +  cos  /3  (sin2  y  -  sin2  /3  ) ^ / sin  y 

2  2  y  2 

=  cos  xp  (1  —  sin  y  sin  0)  2  +  sin  y  sin  0 

2  2  V 

tan  0  =  sin  /3  /  (sin  y  -  sin  /3 )  2  =  sin  0e  /  ( cos  0e  -  sin  y) 


tan  /3  =  sin  y  sin  0/(1  -  sin  y  cos  0  )  =  - ^ ^ - . — 

/  r  .2  .  2  ,*/ 

1 1  —  sin  y  sin  yy  )'2 


tan0e  = 


r  2  2  1/  n  sin  p 

I  [cos  /3  —  (sin  y  —  sin  /3  )'2  ] 

[sin2  /3  +  cos  ^3  (sin2  y  —  sin2  /3  )  ^  ] 

[ ( 1  —  sin2  y  sin"  '*)*  —  s in  y  cos 


=  sin  xp 


[  cos  —  sin2  y  sin2  +  sin  y  sin  ] 
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APPENDIX  V 
FLOW  CHART 
FOR 

SUBROUTINE  ALBEDO  (ALBDO) 


A  E  DC-T  R-65-202 


APPENDIX  V 

FLOW  CHART  FOR  COMPUTER  PROGRAM 


FC-1 


To  Case  Select  in  ($n 


•  a3 
a4 


1  +  sin’1  (Ml\ 

2  \tan  9S  j 

2 tt -  a3 


Eclipse,  No  Albedo  -  Planet 
Radiance  Only 


Subroutine  "ALB  DO"  Arguments: 

0N,  P|S|,  es,  r;  6a 
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Topside,  No  Albedo, 
No  Planet  Radiance 


V  A 

Subroutine  BDI,  with  Arguments  a,  aN,  p,  p,  y,  0.;  Sum,  6a 
Evaluates  the  First  Integral  Between  p,  p 
Adding  Results  to  Sum,  Advancing  a  by  6a 


Subroutine  CXI,  Arguments  a,  aN,  (3^ 
Evaluates  p  on  the  Curve  C.  (a,  p) 
defined  by 

cos  (a  -  aw)  *  -  ctn  p  ctn  Pm 
or 

P  ■  tan-l  r _ li _ : 

[tan  PN  cos  (a  -  aN) 


Is  Test  Split? 


Call 

Subroutine 

BDI 


/ 


Call 

Subroutine 

BDI 


K»2 

a  «  04  +  6 a /2 
a  lim  =  2 7r 


Evaluating 
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To  P-FC-10 
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Subroutine 


AEDC-TR-65-202 
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APPENDIX  VI 
FORTRAN  LISTINGS 
SUBROUTINE  ALBDO 
SUBROUTINE  BDI 
SUBROUTINE  XI 
SUBROUTINE  THET 
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SUBROUTINE  ALBDO ( SUM , BE1  AN , GAM , ALN » THETS* DALS * SGAM » CGAM » STHS , CTHS , 
?  SBN.CBN) _ 


1 

PY  =  3,1415927 

P Y2=  1.5707963 

EP  =  .000001 

2 

SUM  =  0. 

c 

START  THETA  SEARCH 

I F ( THETS  -GAM >400,401, 401 

400 

AL3  =  PY 

AL4  =  PY 

K  =  1 

GO  TO  408 

- 

401 

IFfTHETS  -PY2  >402,402,405  _  __ 

402 

K  =  2 

GO  TO  407 

405 

I F ( THET S  -PY  +GAM  >406,575,575 

406 

K-3 

407 

CAL3  =  -SGAM  *  CTHS/ (CGAM  *  STHS) 

AL3  =  PY2  -  AS  I NF ( CALS ) 

AL4  =  2.*  PY  -  AL3 

c 

408 

START  BETA  SEARCH 

IFfBETAN  -  PY2  +  GAM  >409,409,410 

409 

AL1  =  ALN  +  PY 

AL2  =  ALN  -  PY 

GO  TO  417 

410 

411 

IFfBETAN  -  PY21411, 411*414 

K  =  K+3 

GO  TO  416 

- 

414 

IFfBETAN  -  PY 2  -GAM ) 4 1 5 , 57 5 , 575 

415 

K  =  K+6 

416 

CDAL  =  -CGAM  *  CBN/  (SGAM  *  SBN ) 

DAL  =  PY2  -  AS  I NF ( CDAL  > 

AL1—  ALN  +  DAL 

AL 2  -  =  ALN  -  DAL 

- 

417 

I F ( AL 1  -  2.  *  PY)419,419,418 

418 

AL1  =  AL1  -  2.  *  PY 

419 

I F (  AL2  >420,421,421 

420 

421 


422 


AL2  =  AL2  +  2.*  PY 
GG=K 

GO  TO (422  , 441 .447, 42 5 >452  >520.434*500*540) ,K 

START  CASE  0-0 

ALPHA  =  .5  *  DALS 
BMIN  =0. 


BMAX  =  GAM 
ALIM  =  2.  *  PY 

423  I F (ALPHA  -  ALIM  >424,575,575 

424  CALL  BDI (SUM, ALPHA, DALS, ALN, BMAX  » BMIN, SGAM »CGAM,STHS,CTHS,SBN, CBN ) 
GO  TO  423 

_ START  CASE  0-1  _ 


425 

ALPHA  -  *  5  *  DA 

BMIN  =  0. 

ALIM  =  2 • *  PY 

426 

I F ( ALPHA  -ALIM 
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427 

IF 

AL  1  - 

AL2  >428,428,429 

428 

IF 

( ALPHA 

-AL 1 ) * ( ALPHA  -AL2 

) 1432,431,431 

429 

I F  < 

[ALPHA 

* ( ALPHA  -AL2) )432, 

430,430 

430  IF((ALPHA  -AL1 )*( ALPHA  -2.*PY  ) ) 432 *43 1 ,43 1 

431  BMAX  =  GAM 
GO  TO  433 

432  CALL  XI ( BMAX  ,SBN*CBN* ALN .ALPHA) 

433  CALL  BDI ( SUM  * ALPH A , DALS * ALN , BMAX * BM I N * SGAM * CG AM  *  ST HS , CTHS >SBN , CBN ) 

GO  TO  426  _  ■ _ 1 _ 

START  CASE  0-2 

BMAX  =  GAM 

ALPHA  =  AL2  +  .5  *  DALS 
I F ( AL 1  -AL2)435,436,436 


C 

434 


433  K  =  1 

_ ALIM  =  2«  *  PY _ 

GO  TO  437 

436  K  =  2 
ALIM  =  AL1 

437  I F ( ALPHA  -  AL I M ) 440 , 43 8 , 438 

438  IF( <-1)439*439*575 

439  ALPHA  =  .  5  »■  DALS _ 

GO  TO  436 

440  _ CALL  XI  (  BMIN*_SBN*  CBN*  ALN.  ALPHA)  __ _ _ 

CALL  BDI ( SUM. ALPHA. DALS, ALN, BMAX. BM IN. SGAM *CG AM* SThS* CTHS, SBN*CBN) 

GO  TO  437  _ _ : _  _ _ 

C  ’  START  CASE  1-0 

441  ALPHA  =  .5  *  DALS 
BMIN  _ 

ALIM  =  2e*  PY 

442  I F ( ALPHA  -  AL  IM ) 443  *  575  »  57  5 

443  I F ( ( ALPHA-AL3 )  * ( ALPHA-AL4 ) ) 444 , 445 , 445 

444  CALL  THET  (  BMAX  .STH’SiCTHS  *5GAM  *  ALPHA  ) 

GO  TO  446 

445  BMAX  =  GAM 

446  CALL  BDI ( SUM , ALPHA , DALS , ALN , BMAX , BM I N , SGAM , CG AM , STHS , CTHS ,SBN , CBN ) 
GO  TO  442 

C  START  CASE  2_-0_ 

447  K= 1 

_ ALPHA  =  .5  *  DALS _ 

ALIM  =  AL3 

BMAX  =  GAM  .  _ 

448  I F ( ALPHA  -  AL I M  ) 449 , 450 , 450 

449  CALL  THETIBMIN, STHS, CTHS, SGAM, ALPHA) 

CALL  BDI (SUM, ALPHA, DALS, ALN, BMAX , BMIN, SGAM, CGAM, STHS, CTHS, SBN, CBN) 
_ GO  TO  448 _ 

450  IF(K-1)451»451,575 


451  K  =  2  _ _  _ 

ALPHA  =  AL4  +  .5  *  DALS 

ALIM  =  2.*  PY  _  _  . _ ' 

GO  TO  448 

C _ START  CASE  1-1 

452  BMIN  =  0. 

ALPHA  =  • 5*DALS 
ALIM  =  2 • *P Y 

453  I F ( ALPHA  -AL IM  )  454 , 575 , 575 
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454 
45  5 

456 

457 

458 

459 

460 

461 

462 

463 

464 

465 

466 

467 

468 

469 

470 

471 

472 

473 

474 

475 

476 

477 

478 

479 
48 

48 T  ’ 

482 

483 

484 

485 
603 


601 

604 

486 


487 

602 

488 

489 

490 

491 

C 

500 

501 


I  F ( AL 1  -AL2  )  465  >455  ,455 
I  F ( A  L 1  -AL4)462 ,462 ,456 
.  IF( AL2-AL3 )460  ,460,45  7 

I F (  (ALPHA  “AL3)*(ALPHA  -AL 1 )  )  458 , 602 > 602 
I  F  (  (  ALPHA  -AL2HMALPHA  -AL4  )  )  60  1  >601  *  459 
I F (ALPHA  -AL2 )6C4,603 ,603 

I F (  (ALPHA  — A  L  2 ) * ( ALPHA  -AL 1 )  )  46 1 >602  *  602 
I F ( (ALPHA  -AL3 )* ( ALPHA-AL4) ) 60 1,603, 60 3 
I  F (  (ALPHA  -AL2)*(ALPHA-AL4)  )463>602>6G2 
I F ( (ALPHA  -AL3 ) * ( ALPHA  -A LI) )601 *601*464 
IF (ALPHA  “AL3 ) 603 >603 >604 
IF ( AL 1  -AL3 )476 ,466 ,466 
I  F ( AL 1  -AL4 ) 46  7  ,473,473 
IF ( AL 2  -AL4) 471  ,468 ,468 

I  F (  ( ALPHA  -AL  3 ) * ( ALPHA  “AL2  )  ) 469 > 603  *  603 
IF ((ALPHA  -AL 1 ) * ( ALPHA  -AL4 )  )604, 604*470 
I  F  (  ALPHA  -ALU6C1, 602*602 

IF(  (ALPHA  -AL  3 ) # ( ALPHA  -A L4 )  )  47 2 , 603 > 60 3 
I F (  ( ALPHA  “AL 1 ) * ( ALPHA  -AL2  )  )  601  ,601 , 604 
I F (  (ALPHA  “AL  3 ) * ( A  LPHA  -AL2 )  )  474 > 603  > 603 
I F (  ( A  LPHA  -AL 4 ) * ( ALPHA  - AL  1  )  )  6 03 , 6 03 > 47 5 
IF (ALPHA  -AL4) 601 ,602,602 
IF ( AL 2  -AL3 )483  ,483,477 
IF ( AL2  -AL4)478 ,478,481 

I F (  ( ALPHA  “AL 1 ) * ( ALPHA  -A L4 )  )  479 , 6 03  * 603 
I F (  (ALPHA  “AL 3 ) * ( ALPHA  -AL2 )  ) 6 04 , 480 > 48 C 
I F ( A L PHA  “AL3 ) 6 02*602 *60 1 
IF  ("(ALPHA  “AL1 )  *  (ALPHA  -A  L2  )  )  48  2  >  6  03  *  603 
IF(  (ALPHA  “AL  3 ) * ( AL  PH A  -A L4 )  )  604 , 602 > 602 
I F (  (ALPHA  "AL1 ) * (ALPHA  -AL4 )  ) 48 4 , 603  *  60 3 
I F (  ( ALPHA  -AL2)*(ALPHA  -AL3  )  )  603  *  603  *  48 5 
I  F ( ALPHA  “AL2 ) 602  >602 ,601 

K  =  3 _ 

BMAX  =  GAM 
GO  TO  491 
K  =1 

GO  TO  486 
K  =4 

CALL  THET ( BMAX * STHS , CTHS , SGAM * ALPHA ) 

I  fTk “4 ) 49 1 , 487 *  491 
B 1 =  BMAX 
GO  TO  488 
K  =  2 

CALL  XI (BMAX, SBN, CBN, ALN, ALPHA) 

I F ( K-4 )491 *489,491 
IF (BMAX  -B 1)491  ,491,490 
BMAX  =  B 1 

CALL  BDI ( SUM, ALPHA, DAL S, ALN , BMAX  > BM IN, SGAM, CGAM, STHS, CTHS* SBN ,CBN ) 
GO  TO  453 


ALIM  =  AL 1 _ 

I  F ( AL 1  -AL 2)501*501*502 
K  =1 

ALPHA  =  .5  *DALS 
GO  TO  503 


START  CASE  1-2 
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502 

503 

K  =2 

ALPHA  =AL 2  +  .5*  DALS 

I F ( ALPHA  -ALIM)506>504,504 

5C4 

I  F  €  K  — 

1)505,505,575 

50  5 

ALIM 

=  2.  *  PY 

60  TO 

502 

5C6 

I F (  (ALPHA  -AL 3 ) * ( ALPHA  -AL4 )  ) 508 > 507 : 

507 

BMAX 

=  GAM 

_  _GO_.  TQ 

509 

508  ~ 

CALL 

THET ( BMAX > STHS »CTHS >SGAM> ALPHA) 

509 

CALL 

XI  ( BMIN.SBN.CBN, ALN .ALPHA) 

CALL  801 (SUM, ALPHA  *  DAL S» ALN, BMAX ,  BM IN ,SGAM»CGAM, STHS , CTHS , SBN  ,CBN  ) 
GO  TO  503 


C 

_  52  0 


521 

522 

523 


524 

525 
326 

527 

528 

529 

530 


540 

541 


542 

543 

544 

545 

546 

547 

548 


549 

550 


ALIM  =AL3  ' 

ALPHA  =  .5  *  DALS 
I F ( ALPHA  -ALIM)524»522,522 
IF (<-1)523,523,575 
K  =2 

ALIM  =  2.  *  PY 
ALPHA  AL.4  ’  +"  •5*6ALS 
GO  TO  521 

I F ( AL1  -  AL2  )525,525,526 

I F ( ( ALPHA  -AL 1 ) * ( ALPHA  -AL2 ) ) 52 9 , 5 28 , 528 
I F ( ALPHA  * (ALPHA  “AL 2 ) ) 5 2 9 , 5 2 9 , 52 7 
I F (  ( ALPHAt  J\L1 )*( ALPHA  -2. *P Y ) )  529 > 528, 528 

BMAX  =  GAM  "  . .  . . 

GO  TO  530 

CALL  XI (BMAX, SBN, CBN, ALN, ALPHA) 

CALL  THET (BM IN, STHS,CTHS,SGAM, ALPHA) 

CALL  BDI ( SUM, ALPHA, DALS, ALN, BMAX ,BM I N , 5GAM , CGAM , STHS , CTHS , SBN , CBN ) 
GO  TO  521 

START  CASE  2-2 


BMAX  =  GAM 

I F ( AL 1  -AL2)541,545,545 

K  =1 

ALPHA  =  .5  *DALS 
I  F  (  AL  1  r_AL3J54  2 ,542 ,543 
ALIM  '=  A  L 1  ““ 

GC  TO  553 
ALIM  =  AL3 

I F ( AL 2  -AL4) 544,553,553 
K  =3 

GO  TO  553 _ _ _ 

K  =2 

I F ( AL 1  - AL3  )  5  5 1  ,551,546 
I F ( AL 1—  AL4) 547 ,547 ,549 
I F ( AL 2  -AL3 )548 ,575 ,575 
ALPHA  =  AL2  +.5*  DALS 

A  LJLM_  A  L  3_  _ _ 

GO  TO  553 

I F ( AL 2  ~AL4)550,550,551 
ALPHA  =  AL4  +.5*  DALS 
GO  TO  552 
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551  ALPHA  =  AL?  +  .5  *DALS 

552  ALIM  =  AL  1 

553  I F (ALPHA  -ALIM) 556  ,5  54  ,5  54 

554  I F ( K-2 ) 555  ,575.555 

555  ALIM  =  2.*  PY 
K  =2 

I F  (  AL 2  -AL4)556,557,557 

556  ALPHA  =  AL4  +  . 5*DALS 

_ GO  TO  553 

557  ALPHA  =  AL2  +'.5  *DALS 

GO  TO  553 

558  CALL  THET ( BM I N , STHS , CTHS , SGAM  , ALPHA ) 

IF(K-3)559,561 .559 

559  CALL  X  I  (  B 1  , S9N  , CBN  ,  ALN  ,  AL PHA  ) 

_ I  F( B1  -BMIN ) 561 .561,560 

560  BMIN  =  §1 

561  CALL  BDI ( SUM » ALPHA » DALS , ALN , BMAX »BM I N »SGAM»CGAM »STHS ,CTHS »SBN *CBN  ) 
GO  TO  553 

575  SUM  =  SUM  *  DALS  /  ( SGAM  *4.  *  PY) 

I F ( SUM  -  .00001)576,577,577 

576  SUM  =  L. 

577  BETAN  =  GG 
RETURN 


END 
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1 


3 

4 

5 

6 

7 

8 


9 


10 

11 

12' 

3  3 
14 


2 


SUBROUTINE  BDI  (  SUM  *ALPHA,DALS,ALPHN»  BMAX  »BMIN>SGAM*CGAM>STHS,CTHS  , 


SBETN, CBETN) 

I F (BMIN 

-BMAX 

)1»2,2 

CALPD  = 

COSRF 

(ALPHA  -ALPHN ) 

CALPH  = 

COSRF 

(ALPHA) 

CBMAX  = 

COSRF 

(  BMAX) 

SBMAX  = 

SINRF 

(  BMAX  ) 

CBMIN  = 

COSRF 

(BMIN) 

SBMIN  = 

SINRF 

(  BMIN ) 

ARG1  =  1  •-  ( SBMAX  *SBMAX  )  /  (  SGAM  *  SGAMr) 
ARG2  =  1 • “ ( SBMI N  *SBM I N ) / ( SGAM  *  SGAM) 


I F ( ARG1 ) 3  » 3  >4 
CPMAX  =0. 

G  0  T  0  5  _ _ 

CPMAX  =  SQRTF(ARGl) 

I  F ( ARG2 ) 6  >  6  *7 
CPMIN  =0. 

GO  TO  8 

CPMIN  =  SORTF ( ARG2 ) 

ARG3  =1.  “CPMAX  *CPMAX 
A RG4  = 1 •  -CPMIN  *CPMIN 
I F ( A  RG3 ) 9  *  9  > 1 0 
PM AX  =0  * 

GO  TO  11 

PM AX  =  ASINF ( SQRTF( ARG3 ) ) 

I F ( A  RG4 )12»12»13 

~  "PM  I N  =0  «  ~ . ~ 

GO  TO  14 

PM  I  N  =  AS  I NF ( SQRTF ( ARG4 ) ) 

FE1  =  SBETN*  STHS  *  CALPH  *  CALPD 

FE2  =  CBETN  *  CTHS 

FE3  =  SBETN  *  CTHS . *  CALPD 

FE4  = 'CBETN  *  STHS  *  CALPH 

51  = ( 2 • *  FE1  -.5*(  FE1 _ +  FE2 ) *  CGAM*CGAM  )*  (SGAM*  (CBMAX  *CPMAX 

1  -  CBM  IN  *CPM  I N  )  -  CGAM  *CGAM  *"LOG'fTC  CBMAX+  SGAM*  CPMAX)/(  CBM  I N 

2  +  SGAM  *  CPMIN)))  -  ( FE2  +  FE1)  *5GAM*SGAM*SGAM  *( CBMAX*  CPMAX* 

3  CPMAX  *CPMAX  -  CBM  I N  *  CPMIN  *  CPMIN  *  CPMIN) 

52  =  ( FE4  -FE3 )  *  (.SBMAX*  CPMAX  *CPMAX  *CPMAX  -  SBMI  N  *  CPMIN  *CPMIN 

1  *C  PM  I N  -.5  *  (  SBMAX  *CPMAX  -  SBMIN  *CPMIN  +SGAM  *(PMAX  -PM  IN))) 

2  *SGAM*SGAM*SGAM  _  _ _ _ 

53  = ( FE2  +FE1 ) *  (  SBMAX**4  -  SBMIN**4  ) 

_S4  =  (.5  *FE4+_  1.5  *FE3 )  *  ( BMAX  -BMIN  -SBMA’X  «CBMAX  +SBMIN  *CBMIN 
1  )  +IFE4  -FE3)*(SBMAX*SBMAX  *SBMAX  *CBMAX  -SBMI N*SBMI N*SBM IN*CBM I N ) 

DSUM  =  SI  +S2  +S3  +S4 _ 

SUM  =  SUM  +  DSUM 
ALPHA  =  ALPHA  +  DALS 
RETURN 
END 
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SUBROUTINE  X  I  ( B * SBN , CBN , AN t A L) 

B  =  0. _ 

TB~  =  -CBN  /  (SBN  *  COSRF(AL  -  AN)) 

B  =  ATANF(TB) 

RETURN 

END 

SUBROUTINE  THET ( B  >STHS  * CTHS  * SGAM > AL ) 
F_E_=^CTHS  /  (STHS  *  COSRF(AL)) 

B  =  0. 

B  =  ATANF(FE) 

FE  =  SINRF (B ) / ( COSRF (B )-SGAM) 

FE  =  ATANF(FE) 

B  =  FE  -B 

RETURN _ 

END 
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